This code was executed on the results after running DEEPEst with 5000 iterations on CCP0, CCP1, CCP3, and CCP4 (the 4 studies in the current version of the paper for which DEEP estimates are important; Studies 1-4).
Load data and packages
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 3.5.3
## This is bayesplot version 1.7.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
library(DEEPEst)
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For improved execution time, we recommend calling
## Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
## although this causes Stan to throw an error on a few processors.
load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp1_EstHier8questions.Rdata")
c1est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp1_EstHier8questions.csv")
ggplot(c1est, aes(x = beta)) + geom_density()
ggplot(c1est, aes(x = delta)) + geom_density()
Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP1 (barely). Antonia said rhat should be under 1, which we don’t quite pass.
#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.026951
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.011547
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.145379
#beta and delta (annual discount factor) correlation
cor(c1est$beta,c1est$delta)
## [1] 0.2899706
#beta and daily discount rate correlation
cor(c1est$beta,c1est$r)
## [1] -0.1417037
plot(c1est$beta,c1est$delta)
Antonia does this one participant at a time, so it’s not clear to me what the overall takeaway is. Seems like most pairs plots between beta and delta look okay, though beta and r looks less good (perhaps because r has big positive skew?)
First 12 participants shown below
pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))
pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))
pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))
pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))
pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))
pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))
pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))
pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))
pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))
pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))
pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))
pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))
load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp3_EstHier8questions.Rdata")
c3est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp3_EstHier8questions.csv")
ggplot(c3est, aes(x = beta)) + geom_density()
ggplot(c3est, aes(x = delta)) + geom_density()
Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.
#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.018855
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.014837
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.006831
#beta and delta (annual discount factor) correlation
cor(c3est$beta,c3est$delta)
## [1] 0.2901697
#beta and daily discount rate correlation
cor(c3est$beta,c3est$r)
## [1] -0.1098295
plot(c3est$beta,c3est$delta)
First 12 participants shown below
pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))
pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))
pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))
pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))
pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))
pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))
pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))
pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))
pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))
pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))
pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))
pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))
load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp4_EstHier12questions.Rdata")
c4est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp4_EstHier12questions.csv")
ggplot(c4est, aes(x = beta)) + geom_density()
ggplot(c4est, aes(x = delta)) + geom_density()
Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.
#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.016668
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.006234
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.012964
#beta and delta (annual discount factor) correlation
cor(c4est$beta,c4est$delta)
## [1] 0.2832149
#beta and daily discount rate correlation
cor(c4est$beta,c4est$r)
## [1] -0.1538203
plot(c4est$beta,c4est$delta)
First 12 participants shown below
pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))
pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))
pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))
pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))
pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))
pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))
pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))
pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))
pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))
pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))
pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))
pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))
load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp0_EstHier16questions.Rdata")
c0est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp0_EstHier16questions.csv")
ggplot(c0est, aes(x = beta)) + geom_density()
ggplot(c0est, aes(x = delta)) + geom_density()
Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.
#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.050258
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.039307
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.013776
#beta and delta (annual discount factor) correlation
cor(c0est$beta,c0est$delta)
## [1] 0.3862682
#beta and daily discount rate correlation
cor(c0est$beta,c0est$r)
## [1] -0.21745
plot(c0est$beta,c0est$delta)
First 12 participants shown below
pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))
pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))
pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))
pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))
pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))
pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'
pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))
pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))
pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))
pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))
pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))
pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))